Generic Sandpile Models Have Directed Percolation Exponents 
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We study sandpile models with stochastic toppling rules and having sticky grains so that with 
a non-zero probability no toppling occurs, even if the local height of pile exceeds the threshold 
value. Dissipation is introduced by adding a small probability of particle loss at each toppling. 
Generically, for models with a preferred direction, the avalanche exponents are those of critical 
directed percolation clusters. For undirected models, avalanche exponents are those of directed 
percolation clusters in one higher dimension. 
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In recent years, there has been a lot of interest in the 
study sandpile models as models of real granular matter 
jjj, and also as paradigms of self-organized critical sys- 
tems in general Following the well-known work of 
Bak, Tang and Wiesenfeld many different types of 
sandpile models with different toppling rules have been 
studied B : deterministic and stochastic, with or with- 
out preferred direction, different instability criteria or 
particle distribution rules ||, with fixed energy |tJ etc.. 
However, a clear picture of the factors that determine 
different universality classes of critical behavior is yet to 
emerge M. 

A different paradigm for non-equilibrium critical phe- 
nomena has been directed percolation (DP) which is be- 
lieved to describe active to absorbing state transition in a 
wide class of reaction-diffusion systems |J . The activity 
of avalanches in sandpile models can grow, diffuse or die, 
and any stable configuration is an absorbing state. Thus, 
this should be in the universality class of DP with many 
absorbing states jlCfl . Several models of self-organized 
criticality which show critical exponents related to DP 
have been studied earlier [^T). However, these models do 
not involve any conserved field. The critical exponents of 
known models with conservation of sand are very differ- 
ent from those of DP, and this is presumably due to role 
of the local conservation of sand in the model. In || , a 
model with conservation of particles showing DP expo- 
nents was studied, but this study was mainly numerical. 

In this Letter we study several sandpile models with 
stochastic toppling rules, both directed and undirected. 
The grains are 'sticky' in our models in the sense that 
there is a nonzero probability that any grains arriving at 
a site during the avalanche process just get stuck there. 
We find that the distribution function of avalanches in 
these models has the same power law distribution as that 
of the critical DP clusters. Our theoretical arguments, 
supported by numerical simulations, show that generi- 
cally these models belong to the DP universality class. 
The previously studied deterministic models ||[l2), and 
the stochastic toppling model introduced by Manna |f| 
are unstable to perturbations and flow to the DP fixed 



point. 

The relation of sandpile models to DP was attempted 
earlier in fl4|| using sticky grains. However, in that pa- 
per, no bulk dissipation was present, and sand was dissi- 
pated only at the boundaries. The boundaries break the 
translational invariance of the steady state. The density 
becomes space dependent, and the density profile affects 
the statistics of avalanches. The critical exponents of 
avalanches are not those of critical DP clusters, but ex- 
pressible in terms of DP exponents |l5| ] . Introducing bulk 
dissipation in this paper, we are able to get rid of these 
problems, and the relationship between the sandpile and 
DP problems becomes more transparent. 

For simplicity, we start by defining the directed model 
on a (1 + l)-dimensional square lattice. Various general- 
izations to higher dimensions, undirected case, and other 
toppling rules are straight forward, and briefly discussed 
at the end. The sites on an L x M torus are labelled 
by euclidean coordinates (i,j) with (i + j) even and j 
increasing downwards. At each site there is a non- 

negative integer hi j to be called the height of the pile at 
that site. Initially all hij are zero. The system is driven 
by choosing a site at random and increasing the height 
at that site by one. If one or more particles are added to 
a site at time t (from outside or from other sites), and its 
height becomes greater than 1, then it is said to become 
unstable at time t. 

Any site made unstable at time t relaxes at the 
time (i+1) stochastically: With probability (1— p), it be- 
comes stable without losing any grains. We say that the 
added particle(s) sticks to the existing grains. Otherwise 
(with probability p), the relaxation occurs by toppling 
in which the height at the site decreases by two, and the 
site becomes stable. We introduce bulk dissipation, by as- 
suming that in each toppling, there is a small probability 
5 that both grains from the toppling are lost, other wise 
(with probability 1 — (5), the two grains are transferred 
to the two downwards neighbors (i± l,j + 1). 

Note that there is a nonzero probability that a stable 
site can have arbitrarily large heights. We relax all the 
unstable sites by parallel dynamics. A site made unstable 
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at time t is relaxed in one step at time t+1, independent 
of whether it received one or more grains at the previous 
time step. Once a site has relaxed, it remains stable 
until perturbed again by new grains coming to the site. 
This relaxation process is repeated until all sites become 
stable, and then a new grain is added. 

The model is specified by two parameters p and 5. The 
case p = 1,5 = is exactly soluble, and its critical ex- 
ponents are known in all dimensions JL^ ]. In this case, 
one has to introduce open boundary conditions to ensure 
the existence of a steady state. In two dimensions, the 
probability that adding a particle will cause an avalanche 
of s topplings varies as s~ Ts for large s with t s = 4/3. 
The probability that the duration of avalanche equals T 
varies as T~ Tt with r f = 3/2. The case 5 = 0,p arbitrary 
was studied earlier in [fuf discussed in the introduction. 




No Steady State 



0.2 0.4 0.6 0.8 1 

8 

Fig. 1 The line p = p*(S) separating the steady state 
and no-steady-state regions in the p-5 plane. The inset 
shows a plot of [p c — p*(S)] versus 5. The straight line 
shows the theoretical slope 1/7. 

The probability distribution of different configurations 
in the steady state for general values of p and 5 have an in- 
teresting structure. Define a variable gij — hi j (mod 2). 
We group together different stable configurations {hij} 
corresponding to the same values of {gij}- There are 
2 LM such equivalence classes. On addition of a particle 
to a site, the g variable at that site flips. A toppling leads 
to flipping of g's at the two downward sites. From de- 
tailed balance, it follows that in the steady state, each of 
the 2 LM equivalent classes occurs with equal probability 
& 

If p = 1, with 8 arbitrary, the only allowed height val- 
ues are and 1. In this case, we get a full characterization 
of the steady state. The n-point correlation functions 
satisfy linear equations, and exact solution of can be 
generalized to arbitrary 5. We omit the details here. The 
distribution of avalanche-sizes has an exponential decay 
for non-zero 5. 

As p is decreased below 1, heights greater than 1 ap- 
pear with non-zero probability in the steady state. The 
mean height (hi t j) increases as p is decreased, and there 
exists a critical value p*(S) such that for p < p*(S), the 



height of the pile increases without bound, and there is 
no steady state (Fig. 1). 

The absence of a steady state is obvious along the the 
line p = 0, with 5 arbitrary, as in this case, different sites 
cannot influence each other, and each added particle just 
sticks to the existing pile. A similar decoupling occurs 
along the line 5=1, with p arbitrary. In this case, the 
average particle loss per added grain is 2(1 — co)p, where 
Co is the density of sites with height 0. In the steady 
state this must equal 1. As Co > 0, a steady state can 
exist only for p > 1/2. Thus p*(8 = 1) = 1/2. 

We now derive the equation for the boundary line 
p = p*(8). At the phase boundary, clearly cq — 0. In 
the growth of an avalanche in the steady state of the 
system, any site which receives at least one grain from 
its upward neighbors sends grains to downward neighbors 
with a probability p(l — 5). Thus the probability Prob(s) 
that an added grain will cause an avalanche in which at 
least s sites transfer particles to downward neighbors is 
same as the probability Probp)p(s\p) of a cluster of at 
least s sites in a directed site percolation process with 
concentration of active sites= p. A 

Prob(s) = Prob DP (s \ p = p(l — 5)). (1) 

Note that Eq.(l) also holds in the entire 'no steady state' 
phase, where mean height continues to increase, but 
Prob(s) tends to a limiting distribution for large times. 

By particle conservation, in the steady state, the aver- 
age number of topplings in an avalanche must be equal 
to 1/(25). Let n DP (p) is the sum of the average number 
of occupied and perimeter sites in the cluster correspond- 
ing to a randomly picked site in the DP problem. Now, 
the equation for the phase boundary can be expressed in 
terms of the function n DP (p) as 

n DP ((l-8)p*(6)) = 1/(25). (2) 

For small 5, the average size of clusters is large, and 
p is near the critical probability p c for the directed site 
percolation on this lattice, the function np>p is known to 
vary as (p c —p)~ 7 . Substituting this in Eq.(2) we get 

p*(8) = p c — A5 1 /" 1 + terms of higher order in 5, (3) 

where A is some constant, and 7 is the susceptibility 
exponent of the DP problem. In particular, we have 
p*(8 = 0) = p c - The inset in Fig. 1 shows a log-log plot 
of the numerically determined values of p c — p* (8) plot- 
ted versus 5. We get good agreement with Eq.(3) using 
the existing estimates p c = 0.70548515 and 7 = 2.277730 

Since the average number of topplings in an avalanche 
in the steady state is 1/(25), we will see self- organized 
criticality with long-ranged correlations only in the limit 
5 — > + , and p c < p < 1, (marked with a heavy line in 
Fig. 1). For 5 = + and 1 > p > p c , in the steady state 
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Co > 0. Each site is characterized by two probabilities, 
Pi and P2, which correspond to the probability that a 
site topples when one or two particles are added to it 
respectively. Clearly, pi = (1 — co)p and P2 = p. The 
correlations between heights at different sites are small, 
and if we ignore them |H| ], the evolution of avalanches 
(Fig. 2) in the steady state is as in the Domany-Kinzcl 
model of DP with parameters (pi,p 2 ) pl| . 

DP , , ^ 



is transferred to each of the 2d neighbors of the toppling 
site. Clearly, in this case also, there is no steady state 
for small p. For p = p*(6), the mean height per site di- 
verges. Then, any site which has one or more particles 
added to it, sends particles to its neighbors with prob- 
ability p(l — S). If we look at the space-time history 
of the evolution of the avalanche, we get a directed site- 
percolation cluster on a (d+l)-dimensional body centered 
hypercubical (bch) lattice. The phase boundary in this 
undirected model is also the same as that of the directed 
model on the (d + l)-dimensional bch lattice. 
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Fig. 2 Picture of a typical avalanche for the 2-d di- 
rected model (A) and the time-evolution of the 1-d undi- 
rected model (B) for p = .873 > p c and 5 = .0001 
are compared with the clusters at the critical line of the 
Domany-Kinzel model of DP with p 2 = p. 

Even if some short-range correlations are present, they 
should not modify critical behavior, which is expected to 
be same as in DP. In Fig. 3, we have compared the the 
probability distribution of avalanches with that of DP 
clusters. For s > 1, the latter is expected to satisfy the 
scaling form 

Prob(s) « -A^f[s/ a % (4) 

where A is some amplitudes, r is a critical exponent for 
cluster size distribution in DP and the function f(x) 
tends to a finite constant as x tends to zero, and de- 
creases exponentially with x for large x. We assume that 
s* varies as . Using the constraint (s) ~ we 
get <j) = 1/(2 — t). Using the known numerical estimate 
r = 1.108 in d = 1 + 1 jlffl, we find a very good collapse 
when s T ~ 1 Prob(s) is plotted versus s5^ for two different 
values p — p c and p = 0.873, and two values of 5 = 10~ 3 
and 5 — 10~ 4 . In the inset, the scaling function / is 
compared with that for DP clusters. We get an excellent 
collapse, a strong evidence that the two functions are the 
same, and the correlations in heights in the steady state 
are irrelevant. 

It is straightforward to extend the previous discussion 
to higher dimensions. Thus the avalanche exponents in 
the (d + l)-dimensional model are the same as the expo- 
nents of cluster size distribution in the (d+l)-dimensional 
DP at critical point. The upper critical dimension is 
d = 4. 

Consider now the undirected version of the problem 
on a d-dimensional hypercubical lattice. The rules are 
the same as before, except that on toppling, a particle 
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Fig. 3 Scaling collapse of s T_1 Prob(s) versus s5^ for 
four different combinations of parameters p, S. The inset 
compares the scaling function for the directed and undi- 
rected models with (p, 6) = (0.873, 10 -4 ) with that for 
DP at p = 0.7. 

For p > p*(5), the undirected model differs from the 
directed model in that the height of the pile at a site 
does not change between two topplings. This may give 
rise to possible long term memory effects. However, in 
our model this effect is rather small, as the probability 
of toppling depends on height only if the height is zero. 
Along the line p = p* (5) , the memory affect is strictly ab- 
sent, as the density of sites with zero height goes to zero. 
We find that even for p as large as .873, the avalanches 
are qualitatively similar to DP (Fig. 2), and the dis- 
tribution function is also indistinguishable from that of 
near-critical DP clusters (see inset of Fig. 3). There is 
a crossover from deterministic limit (p = 1, 6 = 0) J3j to 
DP for p ^ 1. 

Why does the conservation of particles not change the 
critical behavior away from DP in our problem? Consider 
a simple DP process of particles of type X, on a lattice 
with nx{r) particles of X at site f. We now attach a 
register ny(r) at each site r, which decreases (increases) 
by 1 each time a particle X is created ( destroyed) at f. 
Then clearly, we have a local conservation of nx + ny- 
Clearly, if the dynamics of the X particle is not affected 
by the book keeping, the process still belongs to the DP 
universality class. In our model, ny is the height of the 
pile. It fluctuates about its mean value, but there is an 
influence of ny on the dynamics of X particles as birth 
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of X particles is not allowed if ny is zero. 

The phenomenological evolutions for the coarse- 
grained density fields nx(r) and ny (r) in the conserva- 
tive limit 8 = + may be written as EC 



d t n x = XJ 2 n x 
d t (n x + n Y ) = 



anx—bn x 



cn x g(n Y ) + n(r,t) (5) 
(6) 



where n(f, t) is a noise term, and a, b 6and c are phe- 
nomenological constants. The crucial term here is the 
coupling term cnxg(ny). Vespignani et al ]20[ ] chose 
g{ny) proportional to ny. In our case, ny has a thresh- 
old and the effect on nx saturates to a finite value even 
as ny increases to infinity. The simplest choice of g{ny) 
to model this is to choose g(ny) — B(ny — h c ), where 9 is 
the step function and h c is the threshold value. In naive 
power counting, this term has the same scaling as the lin- 
ear term in n x . We are not able to treat this analytically. 
But the results of our simulations strongly indicate that 
this perturbation does not change the critical exponents 

ill. 
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Fig. 4 A schematic flow diagram of renormalization 
group flows between different fixed points of sandpilc 
models. 

The DP fixed point is expected to be rather robust 
against perturbations. We have tested several variations 
of toppling rules in simulations. One can make the parti- 
cle transfer process stochastic with each transferred par- 
ticle going to a randomly chosen downward neighbor or 
one can allow multiple topplings at a site with each site 
toppling twice, thrice, etc., with decreasing probability, 
so long as the height is > 1 . With both multiple topplings 
and stochasticity in particle transfer, in the limit of no 
stickiness this becomes Manna model For sticky 

grains, in all these models, we get the DP behavior. The 
schematic renormalization group flows are shown in Fig. 
4. 

To summarize, we have studied several sandpile models 
which show DP-like behaviour. A feature, which is com- 
mon in all these models is 'stickiness', i.e. with a nonzero 
probability a site can remain stable even with a height 
greater than the threshold. This behavior seems to be 
robust against perturbations, and is the generic behavior 
of sandpile models, both directed, and undirected. 

We thank M. Barma, R. Dickman and A. Vespignani 
for their comments on the manuscript. 
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